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Abstract 

We present the results of non linear, hydrodynamic simulations, in three di- 
mensions, of the tidal perturbation of accretion discs in binary systems where 
the orbit is circular and not necessarily coplanar with the disc mid-plane. The 
accretion discs are assumed to be geometrically thin, and of low mass relative to 
the stellar mass so that they are governed by thermal pressure and viscosity, but 
not self-gravity. The parameters that we consider in our models are the ratio 
of the orbital distance to the disc radius, D/R, the binary mass ratio, Mg/Mp, 
the initial inclination angle between the orbit and disc planes, S, and the Mach 
number in the outer parts of the unperturbed disc, Ai. Since we consider non 
self-gravitating discs, these calculations are relevant to protostellar binaries with 
separations below a few hundred astronomical units. 

For binary mass ratios of around unity and D/R in the range 3 to 4, we find 
that the global evolution of the discs is governed primarily by the value of 
A4. For relatively low Mach numbers (i.e. = 10 to 20) we find that the 
discs develop a mildly warped structure, are tidally truncated, and undergo a 
near rigid body precession at a rate which is in close agreement with analytical 
arguments. For higher Mach numbers {Ai ~ 30), the evolution is towards a 
considerably more warped structure, but the disc none the less maintains itself 
as a long-lived, coherent entity. A further increase in Mach number to A4 = 50 
leads to a dramatic disruption of the disc due to differential precession, since 
the sound speed is too low to allow efficient communication between the disc's 
constituent parts. Additionally, it is found that the inclination angle between 
the disc and orbital angular momentum vectors evolves on a longer timescale 
which is probably the viscous evolution timescale of the disc. 
The calculations are relevant to a number of observed astrophysical phenomena, 
including the precession of jets associated with young stars, the high spectral 
index of some T Tauri stars, and the light curves of X-ray binaries such as 
Hercules X-1 which suggest the presence of precessing accretion discs. 



1 



1 Introduction 



Protostellar discs appear to be common around young stars. Furthermore recent 
studies show that almost all young stars associated with low mass star forming 
regions are in multiple systems (Mathieu 1994 and references therein). Typical 
orbital separations are around 30 astronomical units (Leinert et al. 1993) which 
is smaller than the characteristic disc size observed in these systems (Edwards 
et al. 1987). It is therefore expected that circumstellar discs will be subject to 
strong tidal effects due to the influence of binary companions. 

The tidal effect of an orbiting body on a differentially rotating disc has been 
well studied in the context of planetary rings (Goldreich & Trcmaine 1981), plan- 
etary formation, and generally interacting binary stars (cf. Lin & Papaloizou 
1993 and references therein). In these studies, the disc and orbit are usually 
taken to be coplanar (sec Artymowicz & Lubow 1994). However, there are ob- 
servational indications that discs and stellar orbits may not always be coplanar. 
For instance, according to Corporon, Lagrange & Beust (1996), the pre-main 
sequence Herbig star TY CrA is a triple system in which the angle between the 
orbital plane of the third component and that of the central binary has been 
inferred from spectroscopic observations to be ~ 70 degrees. In addition a dusty 
shell is observed around the binary (Bibo, The & Dawanas 1992). Although it 
is located at quite a distance from the central binary (up to about 1000 astro- 
nomical units), this shell is believed to be the remnant of an accretion disc. This 
would indicate that discs and orbits are not necessarily coplanar in pre-main 
sequence systems. In addition, reprocessing of radiation from the central star 
by a warped non coplanar disc has been suggested in order to account for the 
high spectral index of some T Tauri stars (Terquem & Bertout 1993, 1996). A 
dynamical study of the tidal interactions of a non coplanar disc is of interest 
not only in the above contexts, but also in relation to the possible existence of 
precessing discs which may define the axes for observed jets which apparently 
precess (Bally & Devine 1994). 

Linear calculations of the evolution of warped viscous discs were performed 
by Papaloizou & Pringle (1983). The related problem of the propagation of 
bending modes in inviscid discs was considered by Papaloizou & Lin (1995). 
Papaloizou & Terquem (1995) investigated the tidal problem for an inviscid disc 
using linear perturbation theory. Their results suggested that the disc would 
precess approximately as a rigid body if the sound crossing time was smaller 
than the differential precession frequency. They also estimated that the tidally 
induced wave angular momentmii flux arising from bending distortions would 
lead to a small effective viscosity which may in some rather extreme cases be 
significant for protostellar discs. 

In this paper we extend the previous work and consider the full nonlinear 
response of a tidally interacting disc which is not coplanar with the binary orbit. 
We study the problem using a Smoothed Particle Hydrodynamics (SPH) code 
originally developed by Nelson & Papaloizou (1993, 1994). In particular we 
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aim to establish to what extent strongly warped but thin discs may survive in 
close binary systems and how the limitation of the disc size through the effect 
of tidal truncation operates when the disc and binary orbit are not coplanar. 
We find that the tidal truncation effect is only marginally affected by lack of 
coplanarity. Also our model discs were able to survive in a tidally truncated 
condition while warped and undergoing rigid body precession provided that the 
Mach number in the disc was not too large (i.e. if the disc was not too thin). 
An additional result is that the inclination of the disc was found to evolve only 
on a long timcscalc likely to be the viscous timescale, as was indicated by the 
linear calculations of Papaloizou & Terquem (1995). 

The plan of the paper is as follows: In section 2 we give the basic equations. 
In section 3 we describe the orbital configuration and in section 4 we discuss 
the precession of the disc in general terms. In section 5 we review the numerical 
method and in sections 6 and 7 we discuss its application to the disc problems 
considered and estimate the amount of shear viscosity present. Numerical results 
for discs with mid-planes inclined at 45 and 90 degrees to the orbital plane are 
presented for ratio of orbital separation to disc size between 3 and 4 in section 
8. Mass ratios of unity and 0.3 have been considered. Finally, we compare our 
results with some recent observations of relevant astrophysical phenomena in 
section 9. 



2 Basic equations 

The equations of continuity and of motion applicable to a gaseous viscous disc 
may be written 

f + ,V.v = 0, (1) 

^ = --VP-Vt + S„„„ (2) 

at p 

where 

d d „ 
— = h V • V 

dt dt 

denotes the convective derivative, p is the density, v the velocity and P the 
pressure. The gravitational potential is "i, and Smsc is the viscous force per 
unit mass. 

For the work described here, we adopt the polytropic equation of state 

P = Kp^ 

where 
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c2 = if7p^-i 

gives the usual associated sound speed, c^. Here we take 7 = 5/3, and K 
is the polytropic constant. This corresponds to adopting a fluid that remains 
iscntropic throughout even though viscous dissipation may occur. This means 
that an efficient coohng mechanism is assumed. 

3 Orbital Configuration 

We consider a binary system in which the primary has a mass Mp and the 
secondary has a mass Mg. The binary orbit is circular with separation D. The 
orbital angular velocity is uj. Wc suppose that a disc orbits about the primary 
such that at time t = it has a well defined mid-plane. We adopt a non rotating 
Cartesian coordinate system (x, y, z) centred on the primary star and we denote 
the unit vectors in each of the coordinate directions by i, j and k respectively. 
The z axis is chosen to be normal to the initial disc mid-plane. We shall also 
use the associated cylindrical polar coordinates (r, Lp, z). 

We take the orbit of the secondary to be in a plane which has an initial 
inclination angle 5 with respect to the (x, y) plane. For a disc of negligible mass, 
the plane of the orbit is invariable and does not precess. We denote the position 
vector of the secondary star by D with D = |D|. Adopting an orientation of 
coordinates and an origin of time such that the line of nodes coincides with, and 
the secondary is on, the x axis at f = 0, the vector D is given as a function of 
time by 

D = £) cos Loti + D sin wt cos 5 i + Dsmujt sin 5 k. (3) 
The total gravitational potential ^' at a point with position vector r is given 

by 

GMp GMs GMsY-T> 
^ = P '- 1 

I r I I r - D I D3 

where G is the gravitational constant. The first dominant term is due to the 
primary, while the remaining perturbing terms are due to the secondary. The 
last indirect term accounts for the acceleration of the origin of the coordinate 
system. We note that a disc perturbed by a secondary on an inclined orbit be- 
comes tilted, precesses and so does not maintain a fixed plane. Our calculations 
presented below are referred to the Cartesian system defined above through the 
initial disc mid-plane. However, we shall also use a system defined relative to 
the fixed orbital plane for which the '.r axis' is as in the previous system and 
the axis' is normal to the orbital plane. If the disc were a rigid body its angu- 
lar momentum vector would precess uniformly about this normal, as indicated 
below. 
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4 Precession Frequency 



In order to aid discussion of our numerical results, we consider in general terms 
the response of the disc to the component of the time averaged perturbing po- 
tential, arising from the secondary, that is odd in z. This is the contribution 
responsible for inducing disc precession. We recall that to lowest order in 1/D 
this is (see Papaloizou & Terquem 1995): 

-r 3GM, 

Ws = — - rz sm 2d smyj. (4) 

We take the unperturbed disc to be in a state of Keplerian rotation with lo- 
cal angular velocity f2. Assuming that the basic equations can be linearized, 
the Lagrangian displacement vector ^ will satisfy an equation of the form (see 
Lynden-Bell & Ostriker 1967) 

C(^) = -V*,. (5) 

Here C is a linear operator, which needs to be inverted to give the response. 
When a barotropic equation of state applies, and the boundaries are free, C is 
self-adjoint with weight p. This means that for two general displacement vectors 
^(r) and J7(r) we have 



JV Vjv 



(6) 



where * denotes complex conjugate and the integral is taken over the disc volume 
V. 

A general time averaged perturbing potential may be Fourier analyzed in 
the direction. We here consider the component with azimuthal mode number 
m = 1. Because of the spherical symmetry of the unperturbed primary potential, 
the unperturbed system is invariant under applying a rigid tilt to the disc. This 
corresponds to the existence of rigid tilt mode solutions to (^) when there is no 
forcing ^ = 0). For a rotation about the x axis, the rigid tilt mode is of the 
form 



^ = 1^ (r X c^) sin — c^z cos (7) 

where (p is the unit vector along the ip direction. For time averaged forcing 
potentials cx sin (/s, the existence of the solution (|^) results in an integrability 
condition for (||). When C is self-adjoint, this is 

/ |-r-V*spdrEE / i- (rxV*,)pdr = 0. (8) 
JV Jv 

The above condition is equivalent to the requirement that the x component of 
the external torque vanishes. This will clearly not be satisfied in the problem 
we consider. 
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To deal with this one may suppose that the disc angular momentum vector 
precesses about the orbital angular momentum vector with angular velocity uSp 
(Papaloizou & Terquem 1995). This in turn is equivalent to supposing our 
coordinate system rotates with angular velocity ujp about the orbital rotation 
axis. Treating the Coriolis force by perturbation theory produces an additional 
term on the right hand side of (||) equal to —2rVLu)pXip. Using the modified 
force in formulating (^) gives the integrability condition as 

ujpSmS / r'^ilpd.T ^ / i • (r X Vl's) pdr, (9) 
Jv JV 

with ujp — \iji3p\. Equation gives a precession frequency for the disc that would 
apply if it were a rigid body. However, approximate rigid body precession is 
only expected to occur if the disc is able to communicate with itself, either 
through wave propagation or viscous diffusion, on a timescale less than the 
inverse precession frequency (see for example Papaloizou & Terquem 1995 and 
below). Otherwise, a thin disc configuration may be destroyed by strong warping 
and differential precession. 



5 Numerical Method 

The set of basic equations (|l]) and (||) are solved numerically using an SPH 
code (Lucy 1977, Gingold & Monaghan 1977), developed by Nelson & Pa- 
paloizou (1994), which uses a conservative formulation of the method that em- 
ploys variable smoothing lengths. A suite of test calculations illustrating the 
accurate energy conservation obtained with this method is described by Nel- 
son & Papaloizou (1994), and additional tests and calculations are presented in 
Nelson (1994). 

5.1 Hydrodynamical Equations 
5.1.1 Smoothed Functions 

In SPH one works with smoothed quantities defined by the expression 

(/(r)>- / f{v')W{\v-v'\,h)d^v' (10) 



where /(r) is an arbitrary function, H^ is a smoothing kernel, and h is the 
smoothing length. Because the function /(r) is defined only at particle positions, 
a Monte-Carlo representation of equation (|l^) is used in the form 



N 

rrii 



/(^^OW^dr.-r.h/i) (11) 

3 = 1 



6 



where rrij is the mass of particle j, pj is the density at the location of particle 
j, and the sum is over all N particles. If the function /(r) is chosen to be the 
density, p{v), then we obtain 

TV 

p(r,)=^™,l^(|r,-r,|,/i) (12) 
i=i 

where the angled brackets have been dropped for the sake of brevity, with the 
understanding that smoothed functions are used throughout in SPH. 

When variable smoothing lengths are employed each particle, i, has an as- 
sociated smoothing length hi. In order to conserve momentum the kernel is 
symmetrised with respect to particle pairs using the procedure of Hernquist & 
Katz (1989) which leads to the replacement: 

W^, ^\[W (|r,; - r,| , h{) + W (|r, - r,| , h,)] (13) 

The smoothing function used is the commonly employed spline kernel ini- 
tially suggested by Monaghan & Lattanzio (1985). 
Then the smoothed density (|l^) becomes 

N 1 

5.1.2 Equations of Motion 

The equation of motion for particle i is 

^7p,i+^G.t + ^visc,i- (15) 

at 

Here we have denoted the forces on particle i due to pressure, gravity and 
viscosity as ¥ps^'PG,i, and S„isc,i respectively. 

5.1.3 The pressure force 

The pressure force is given by the expression 



N 

N N 

k=i j=i 









'dW{r,,,h,) 






5 


d\r,j\ 



dW{r,,,hj) 



hiConst l^'yl 9\rij\ 



fpk^ 


\ 


'dW{vkj,hk) dhk 


\pIj 




dhk dvi 



dvi 



hj const I *i I 

(16) 



where we adopt the notation Tkj = r^ — r^ . A full derivation of this form of the 
pressure force equation is given in Nelson & Papaloizou (1994). The first part 
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of the above pressure force, that does not involve derivatives of the smoothing 
length, is of the usual symmetrised standard form (Hernquist & Katz 1989). In 
order to evaluate the second part, we need to specify the way the smoothing 
length depends on the particle coordinates. 

In order that the system conserve linear and angular momentum as well as en- 
ergy, we take the h^, to be functions only of the absolute distances between parti- 
cles. Here, we shall adopt the prescription given in Nelson & Papaloizou (1994), 
namely we first find the Mtol nearest neighbours and then calculate hi from 
this list of nearest neighbours according to: 

n— 1 

where the summation is over particle i's Nfar most distant nearest neighbours. 
For the calculations presented here, we have used Ntol = 45 and Nfar = 6. 
Then the value of 2hi is equal to the mean distance to particle i's six most 
distant nearest neighbours. 



5.1.4 The addition of viscosity 

In order to stabilize the calculations in the presence of shocks, we add an artifi- 
cial viscous pressure term, 11^ (Monaghan & Gingold 1983), in only the leading 



term of equation (16), such that 



P, P, P, Pi 

9 I 9 ^ 9 ^ 91^ AljT . 
P- P- P- 



With this addition 



Fp^i ^ Fp^ S-uisc^i- 

The artificial viscous pressure given by Gingold & Monaghan (1983) is 

Ily = ^ (-afiijCij + (3fijj) (18) 

Pij 

where the notation Aij = ^(Ai+Aj) has been used, Cy is the mean sound speed 
evaluated at particles i and j, and 

I Ji \ ^2 ifvjj .rij <0 
= ^ (19) 

otherwise. 

Here, the notation v^j = (vi — Vj) has been used, and 77^ = 0.01/ifj- prevents 
the denominator from vanishing. For all calculations reported here, we adopted 
a = 0.5 and /? = (cf. Artymowicz & Lubow 1994). 
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5.1.5 Time stepping 



The particle equations of motion were integrated using tlie standard leap frog 
method. The time step size was determined using the criterion presented in 
Nelson & Papaloizou (1994), with the value of the Courant number being set to 
Q = 0.3. 

6 Initial set up 

A disc containing 17500 identical particles was initially set in orbital motion 
about the primary alone. This was assumed to act through a softened gravita- 
tional potential such that 

where h is the softening length. With this prescription the disc could be extended 
to r = as was the case for the calculations presented here. We adopt units such 
that the primary mass Mp = 1, the gravitational constant G = 1, and the outer 
disc radius R = 1. In these units the adopted softening parameter b = 0.2 and 
the time unit is n{R)~^. The particle positions were initialised by dividing the 
disc into 50 annuli of equal width and placing in each one the number of equal 
mass particles required to maintain a constant E distribution over the whole 
disc, where S is the surface density averaged over an annulus. The particles 
were then separated into seven equally spaced planes, giving a total uniform 
semi-thickness H. 

We note that in standard axisymmetric accretion disc theory, vertical hydro- 
static equilibrium at a general radius gives Cs(r) fl{r)H{r). In the outer parts 
of the disc, this is equivalent to = {H/R)~^, where we define the Mach 
number M as: 

M = v^/cs, 

where the azimuthal velocity v^, and Cg arc evaluated in the outer regions of 
the disc just interior to the disc edge. Hence we say that the Mach number M. 
parameterises the hydrostatic state of the disc. The value of H was estimated 
so as to give vertical hydrostatic equilibrium in the neighbourhood of the outer 
disc boundary for the required Mach number. 

After determining H, the particle positions in each plane section for each annulus 

were generated psciido-randomly by computer. The particle velocities were set 
to be initially azimuthal and equal to rfl, the circular velocity, such that, in our 
units 



(r2 +62)3/4 • 
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Finally the required mean Mach number in the outer regions of the disc, M , was 
made consistent with the equation of state by suitably adjusting the polytropic 
constant K, since at any location = (r'^Vi^)/{K^p''~^). 
For all of the calculations described here the smoothing length was less than 
the disc semi-thickness so that the latter was determined by genuine pressure 
effects. Only for the largest value considered {M = 50) were the smoothing 
length and semi-thickness comparable. 

The disc models were relaxed, before introducing the secondary, for approxi- 
mately two orbital periods at the outer boundary. During this period the outside 
edge of the disc expanded by about 20 percent to 1.2. This effect was mostly 
due to a pressure-driven expansion into the surrounding vacuum, near the outer 
edge, and a small viscous expansion. After relaxation, the Mach number was 
found to be approximately constant through the disc and the expansion at the 
outside was found to proceed through the action of viscosity alone. 

7 Viscosity Calibration 

Although in a formal sense the artificial viscosity specified by (|l8|) is a bulk 
viscosity designed to handle shocks, its practical implementation induces a 
shear viscosity v = C^aCgh, with being a constant usually in the range 
0.1 to 1 (Artymowicz & Lubow 1994). This leads to angular momentum trans- 
port and the standard viscous evolution of an accretion disc (Lynden-Bell & 
Pringle 1974). In studies of the type presented here where discs undergo angu- 
lar momentum exchange through tidal interaction with orbiting secondaries, it 
is important to know the rate at which viscosity transports angular momentum 
so this may be compared with other effects. In the calculations presented here it 
is evident that in many cases, disc expansion arising from outward transport of 
angular momentum is halted by tidal truncation. This effect, well known in the 
coplanar case (Lin & Papaloizou 1993), is seen here when the disc and binary 
angular momenta are not aligned. 

In order to estimate the magnitude of the shear viscosity operating in our disc 
models we studied some freely expanding disc models for up to fifty time units. 
For our variable smoothing length prescription we expect that h cx p^^l'^ . Then, 
as for our equation of state Cs oc p^^^ ^ we expect v to be constant for our disc 
models. To specify the magnitude of v we write v = assC^{R)/^{R), where ass 
corresponds to the well known Shakura & Sunyaev (1973) a parameterization. 
However, it is applicable only at the outside edge of the disc. 

The standard diffusion equation for the evolution of a viscous disc with v a 
constant is (Lynden-Bell & Pringle 1974) 

f^^^{[<^'")'r' I PC) 

where a prime denotes differentiation of a function with respect to r. 
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We found that the evolution of a freely expanding SPH disc, with A4 — 10, 
followed (^0|) to high accuracy for < r < 1 and < i < 50. This was found 
to be the case immediately after set up with O evaluated using the softened 
primary potential, when ass = 0.03 near the outer boundary of the disc. This 
corresponds to a global viscous evolution timescale of 4000 time units. Although 
this is a longer time than any of our models were evolved for, the time required 
to establish tidal truncation after a radial expansion of 20 percent is significantly 
shorter, being about 160 units for M = 20. 

We also comment that a value of ass = 0.03 is about two orders of magni- 
tude larger than that expected to be associated with tidally induced inwardly 
propagating waves (Spruit 1987, Papaloizou & Terquem 1995). Accordingly, it 
is expected that tidal truncation will instead occur through strong nonlinear 
dissipation near the disc's outer edge (Savonije, Papaloizou & Lin 1994). The 
large viscosity of the disc models considered here will damp inwardly propagat- 
ing waves before they can propagate very far. 

8 Numerical Results 

We have considered the evolution of the disc models set up according to the 
procedure outlined above. Details of the models are given in Table 1. In the 
absence of the secondary star the models were characterized only by the mean 
Mach number in the outer regions of the disc, Ai . After the relaxation period 
of about two rotation periods at the outer edge of the disc, the time was reset to 
zero and the secondary was introduced in an inclined orbit, moving in a direct 
sense, crossing the x axis at t = 0. The parameters of the binary orbit are given 
m Table 0. In some cases the full secondary mass was included immediately 
corresponding to a sudden start. However, for strong initial tidal interactions 
such as those that occur when D/R = 3, Ms/Mp = 1, this can result in disrup- 
tion of the outer edge of the disc with a small number of particles being ejected 
from the disc. It was found that this could be avoided by using a 'slow start' in 
which Ms was built up according to 

Ms{t) = Ms (1 - Aexp(-ai)) , 

where we used A = exp(— 0.004), and a — 0.4. Consequently, the total compan- 
ion mass Ms was present in the calculations after a time of w 10. 

Models were continued for the run times indicated in Table 0. Most binary 
systems had D/R — 3, but cases with D/R ~ 3.6 and 4 were also considered. 
The general finding was that a disc with an initial angular momentum vector 
inclined to that of the binary system tended to precess approximately as a rigid 
body, with a noticeable but small warp if M was not too large. In such cases only 
small changes in the inclination angle between the angular momentum vectors 
were found over the run time. This is consistent with the expectation from 
Papaloizou & Terquem (1995) that the timescale for evolution of the inclination 
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in such cases should be comparable to the viscous evolution timescale of the 
disc, assuming outward disc expansion is prevented by tidal interaction. 



Table 1: Model discs 



Model Ms/Mp D/R 5 M Start Runtime 



1 


1 


3 


7r/4 


20 


slow 


310.0 


2 


1 


3 


7r/4 


25 


slow 


217.7 


3 


1 


4 


7r/4 


20 


sudden 


330.0 


4 


1 


3.6 


7r/4 


10 


sudden 


152.0 


5 


1 


3.6 


7r/2 


10 


sudden 


156.0 


6 


0.3 


3.6 


7r/4 


10 


sudden 


84.0 


7 


1 


3 


7r/2 


20 


slow 


328.5 


8 


1 


4 


7r/2 


20 


slow 


228.0 


9 


1 


3 


7r/4 


30 


slow 


397.9 


10 


1 


3 


7r/4 


50 


slow 


50.0 


11 


1 


3.6 


0.0 


10 


sudden 


150.8 



The calculations presented here use a coordinate system which is based on 
the initial disc mid-plane. However, as the disc precesses, the mid-plane changes 
location with time. It is then more convenient to use a coordinate system 
(x, Uo, Zo) based on the fixed orbital plane, the Zo axis coinciding with the orbital 
rotation axis. We shall refer to these as 'orbital plane coordinates'. We locate 
the inclination angle i (equal to (5 at i = 0) between the disc and binary orbit 
angular momentum vectors through 

Jd • Jo 

Here, Jo is the orbital angular momentum. The disc angular momentum is 
Jd = Jj, where the sum is over all disc particle angular momenta Jj. 
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A precession angle (3p, measured in the orbital plane can be defined through 

(J13XJ0) • u 



cos/3p 



|JdxJo||u| 



where u may be taken to be any fixed reference vector in the orbital plane. We 
take this to point along the yo axis such that initially (3p — tt/2 in all cases. For 
retrograde precession of Jd about Jo the angle (3p should initially decrease as 
is found in practice. The period of revolution for rigid body precession of a disc 
is 2-K/bjp with, from equation (||): 

For a disc with constant E and radius R, this gives 



However, near rigid body precession could only be expected if there is good com- 
munication, either through wave propagation or viscous processes, throughout 
the disc. As long as ass < which is at least marginally satisfied in all mod- 

els except model 10, sound waves may propagate (Papaloizou & Pringle 1983). 
The ability of sound to propagate throughout the disc during a precession time 
implies approximately that 

f > ^- 

When Ms/Mp = 1, D/R = 3 and 5 = 7r/4, equation (||) gives ujp/n{R) = 0.012 
corresponding to a precession period of 512 in our time units. Our results were 
consistent with the condition ( |23|) to within a factor of two in that models with 
M < 25,6 = 7r/4 and D/R > 3 showed modest warps and approximate rigid 
body precession. Model 9 with = 30 showed severe warping and a more 
complex precessional behaviour while model 10 with = 50 appeared to be 
disrupted by differential precession. 

We now describe in more detail the results from some of our runs. We 
present particle projection plots for each of the Cartesian planes using obital 
plane coordinates. In all such figures a fourth 'sectional plot' is also included in 
which only particles such that —0.05 < yo < 0.05 are plotted. 



8.1 Models 1 & 3 

A projection plot for model 1 is shown in Fig. (|l|) near t = when the disc is 
unperturbed. Note that the disc appears as edge on and inclined at 45 degrees 
in the (i/o, Zo) plane. The time dependence of the angles t and Pp is plotted 
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in Fig. (^. It may be seen that there is little change in l during the whole 
run. On the other hand (3p decreases approximately linearly, corresponding to 
uniform precession (note that /3p is plotted as positive rather than negative 
for the latter section of this plot). The inferred precession period is around 
600 units, in reasonable agreement with the value of 512 units obtained from 
equation (22). Fig. (||) shows a projection plot at time 156.6 after the disc 
has precessed through about 90 degrees in a retrograde sense. This amount 
of precession is evident because at this time the projection in the (a;, z^) plane 
shows the disc to be almost edge on, and inclined at about 45 degrees. Note 
that the disc remains thin. Fig. (||) shows another projection aX t = 228.8. This 
is taken at a viewing angle such that the disc is almost edge on. Particles from 
the disc at time t = are also projected. It is apparent that at this time the 
disc plane is approximately orthogonal to the original one. The later disc is 
also clearly of smaller size demonstrating the effect of tidal truncation. We note 
that the same configuration, viewed from an angle at which the original disc 
appears edge on, indicates that the apparent smaller size of the evolved disc is 
not the result of projecting the short axis of an elliptical disc. The mean radius 
of the evolved disc is indeed smaller than the original disc radius due to tidal 
truncation. 

Fig. (|^) shows a plot similar to Fig. (|3|) at t = 297.8, near the end of the 
run when the disc has precessed through about 180 degrees. This amount of 
precession is demonstrated by the fact that the disc appears almost edge on in 
the (2/0, Zo) plane just as it did at time t = 0. However, its plane is inclined at 
about 90 degrees to the original disc plane. At this stage our results indicate 
that the disc has attained a quasi-steady configuration as viewed in a frame that 
precesses uniformly about the orbital rotation axis. In this run, as in others, 
the disc develops a warped structure that initially grows in magnitude but then 
levels off. The sectional plot in Fig. (^) indicates that the disc has developed a 
modest warp in this case. 

For comparison with model 1, Fig. (^) shows a corresponding projection for 
model 3 with D/R — A at t = 330.0. As expected the precession rate is about 
half as fast in this case, but the sectional plot indicates that the warp is not 
less pronounced than in model 1. This may be related to the fact that because 
of the more distant companion, the disc is able to expand to a larger size. We 
also note that the linear calculations of Papaloizou & Terquem (1995) indicated 
that the low frequency components of the perturbation do not diminish in their 
effectiveness until D/R ~ 8. 



8.2 Models 9, 2, & 10 

In contrast to models 1 and 4, in which the Mach numbers are Ai = 20 and 
= 10 respectively, the behaviour of model 9 with M = 30 but similar orbital 
parameters is considerably more complex. This disc develops a strong warp such 
that the inner and outer parts of the disc try to separate. A projection plot is 
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shown at t = 397.9 in Fig. (^. The inner part of the disc seems to occupy a 
different plane from the outer part. The outer part was found to precess like a 
rigid body at the expected rate, and it tended to drag the inner section behind 
it. This is indicated in Fig. (||) where we plot the angle f3p calculated using the 
outer disc section with r > 0.5 only and also the same angle calculated using 
only the inner section with r < 0.3. These do not correspond exactly but the 
angle associated with the inner segment progresses at a variable rate indicating 
coupling to the outer section. The angle t associated with the two sections is 
also plotted. This becomes significantly smaller for the inner segment consistent 
with the existence of a large amplitude warp. We note that the relatively larger 
inclination associated with the outer segment enables a closer matching of the 
precession frequencies associated with the two segments and aids coupling, due 
to the presence of large pressure gradients induced by the strong warping. Each 
segment of this model remained thin throughout the run. 

For comparison with model 9, Fig. (|^) shows the evolution of Pp and t for 
model 2, which is the same as model 9 but with Ai = 25. This is an intermediate 
case between model 9 and model 1. The precession is slightly differential at the 
beginning of the run, but the inner part of the disc remains coupled to the outer 
part, and the precession becomes uniform after about 150 time units. Since the 
inner part adjusts its precession frequency to the outer part, it changes slightly 
its inclination. 

In the case of model 10, which had a Mach number of Al = 50, severe disruption 
of the disc due to differential precession was found to occur. The thickness of the 
disc was such that efficient communication between different parts of the disc 
was unable to occur, and consequently a well defined, coherent warped structure 
could not be maintained. Instead, a decidedly non disc-like morphology was 
obtained by the time the calculation had been evolved for t = 50 units. 

8.3 Models 4 & 11 

We also studied relatively thick disc models with M = 10. For example model 4 
had D/R = 3.6, and S = 7r/4. These models also precessed like rigid bodies with 
barely discernable warps. We plot the behaviour of i and /3p (calculated for the 
whole disc) in Fig. (pi]). In this case the precession period is estimated to be 873 
units, again consistent with the value of 884 units predicted by equation (p2|). 
We note that the agreement between theory and the numerical calculations is 
better in the = 10 than for the higher Mach number cases calculation because 
the larger sound speed leads to improved communication throughout the disc, 
allowing it to precess more as a rigid body. Here t shows little evolution over the 
duration of the run. A projection plot is shown at t — 105.1 in Fig. (^0|). This 
figure illustrates truncation of the disc at i? = 1 as a freely expanding disc was 
found to expand significantly further. For comparison we ran the same model 
but with 5 — (coplanar). A projection plot for this model 11 at i = 147.3 is 
shown in Fig. (|lj). It was found that the discs in models 4 and 11 are of about 
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the same size and thickness. Thus the tidal truncation radius for (5 = 7r/4 is 
similar but slightly larger than in the coplanar case. 

8.4 Model 6 

We also ran model 6 which had the same parameters as model 4, but with 
a mass ratio of Ms/Mp = 0.3. The lower secondary mass results in a bigger 
disc, and a precession rate about three times slower, as one would predict from 
equation (p^). A projection plot at i = 84 is given in Fig. (|l3|), in which there 
is little sign of a warp. 

8.5 Models 5, 7 & 8 

In addition to the cases described above we also ran models with 5 = 7r/2. In 
this case the disc plane remains in a stable position apart from some oscilla- 
tory bending and rotation. This is again consistent with inclination evolution 
occuring on the global viscous timescale. These discs also remained thin and 
exhibited tidal truncation, though the discs were able to expand to a somewhat 
larger size in these cases. These features are indicated in the projection plots 
for model 7 at t = 328.5 and model 8 at i = 228.0 shown in Fig. (|l|) and @ 
respectively. Model 8 had D/R = 4 and thus expanded to a larger radius until 
tidal effects became significant. 

9 Discussion 

In this paper we have considered the nonlinear response of an accretion disc to 
the tide of a binary companion when the disc mid-plane is not necessarily copla- 
nar with the plane of the binary orbit. For our constant viscosity SPH models we 
found that tidal truncation operates effectively when the disc and binary orbit 
are not coplanar, being only marginally affected by the lack of coplanarity. Our 
results indicate that modestly warped and thin discs undergoing near rigid body 
precession may survive in close binary systems, with the integrity of their local 
vertical structure maintained. However, extremely thin discs may be severely 
disrupted by differential precession and their survival depends only on their hy- 
drostatic state (i.e. the Mach number, M., of the unperturbed outer disc). For 
the binary mass ratios and separations considered here, the crossover between 
obtaining a warped, but coherent disc structure, and disc disruption, appears to 
occur for values of the Mach number M > 30. We also found that the inclina- 
tion evolved on a long timescale, likely to be the viscous timescale, as indicated 
by the linear calculations of Papaloizou & Terquem (1995). Below, we discuss 
some of the potential applications of our models to a number of astrophysical 
phenomena of current interest. 
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9.1 Precessing jets in star forming regions 

Amongst the currently popular models formulated to explain the generation of 
jets in young stellar objects, there are a class of models that rely on the emission 
of a wind from the disc surface which is then accelerated and coUimatcd by the 
action of a magnetic field (see for example Konigl & Ruden 1993 and references 
therein). It is reasonable to assume that the wind will largely be emitted in 
a direction normal to the disc surface, so that a processing disc may lead to 
the excitation of a precessing jet. The precession periods obtained from our 
calculations with a mass ratio of 1 vary between 500 and 900 in our units. 
When scaled for a disc of radius 50 AU, surrounding a star of IMq, the unit 
of time is rt{R)~^ ~ 56 yr, leading to a precession period of between 3.10'* and 
5.10V- 

Bally & Devine (1994) suggest that the jet which seems to be excited by the 
young stellar object HH34* in the L1641 molecular cloud in Orion processes 
with a period of approximately 10^ yr. This period is consistent with the source 
being a binary with parameters similar to those we have used in our simulations 
with a separation on the order of a few hundred astronomical units. We note 
that such binaries cannot be resolved in Orion (at a distance of 470 pc) with 
current ground based technology, but could be resolved with the Hubble Space 
Telescope. 

9.2 Spectral energy distributions of young stellar objects 

Some of the results presented in this paper demonstrate how a warped disc can 
present a large surface area for intercepting the primary star's radiation. The 
effect that the consequent reprocessing of the stellar radiation field can have on 
the emitted spectral energy distribution has been investigated by Terquem & 
Bertout (1993, 1996). They find that the additional shrouding of the central 
star by a strongly warped disc may account for the high spectral index of some 
T Tauri stars or even for the spectral energy distribution of some class O/I 
sources. Until now it was not clear to what extent a strongly warped disc 
configuration could be created and maintained over long time periods. However, 
the behaviour of Model 9 demonstrates that such a scenario could be physically 
realisable. 

We note that the gradual warping of an accreation disc will lead to it being 
progressively heated by the radiation from the central source, so that the Mach 
number will then be a function of time. Since the Mach number is the deter- 
mining parameter in the global evolution of the disc, we anticipate that the 
inclusion of this effect may lead to changes in the evolution of the discs, par- 
ticularly in the case of model 9, where strong warping was observed to occur. 
Future calculations will be presented to address this and related issues. 
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9.3 Light curves of X— ray binaries 



There is evidence from the hght curves of X-ray binaries such as Hercules X-1 
and SS433, that their associated accretion discs may be in a state of precession 
in the tidal field of the binary companion. We note that these systems are more 
complicated than those modelled here because of the process of mass transfer 
occuring in them. In addition, numerical investigations on the accretion disc 
coronae (Ko & Kallman 1994) suggest disc mach numbers which our results 
imply are marginal for survival under differential precession in a tidal field. 
However, it is still of interest to compare the precession rates that are inferred 
observationally to those that we would predict from our models. We assume 
that the discs around these objects are tilted because of some internal process, 
for example the torque exerted by a coronal wind (Horn & Meyer 1994), and 
that the secondary component induces a precession of the disc. 

From equation (p^), the ratio of the precession frequency to the disc outer 
edge frequency is proportional to (R/D)^ cos (S) /q, where q — Mp/Mg. To 
obtain a general scaling we assume that the disc radius is a fixed fraction of 
the mean Roche lobe radius. The ratio of the Roche lobe radius Rl to the 
separation is a function only of q. We then find the approximate scaling for 



where Porb is the orbital period of the binary. We calculate the factor of propor- 
tionality by using our numerical results. We find that Pprec/ Port varies between 
10 and 20 as q increases from to 1 for small inclinations 5. Observations give 
a ratio close to 20.6 for Her X-1 {Porb — 1-7 days and Pprec = 35 days, Petter- 
son 1975, Gerend & Boynton 1976) and about 12.6 for SS433 (Porb = 13 days 
and Pprec = 164 days, Margon 1984 and references therein). These values are 
consistent with the precession being induced by the tidal field of the secondary. 
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Figure 1: Projection plots in orbital plane coordinates for model 1 at time 
t ~ 0. The projections are in the {x,yo) plane (top left), {x,Zo) plane (top 
right), (ucZo) plane (bottom left) and {x,Zo) plane (bottom right). 

Figure 2: The precession angle f3p (dashed line) and inclination angle l (solid 
line) for model 1. 

Figure 3: Projection plots in orbital plane coordinates for model 1 at time 
t = 156.6. 

Figure 4: Projection plot in the initial disc plane {x, y) for model 1 at time t = Q 
and t = 222.8. The disc at t = 222.8 is seen almost edge-on. 

Figure 5: Projection plots in orbital plane coordinates for model 1 at time 
t = 297.8. 

Figure 6: Projection plots in orbital plane coordinates for model 3 at time 
t = 330.0. 

Figure 7: Projection plots in orbital plane coordinates for model 9 at time 
t = 397.9. 

Figure 8: The precession angle /3p (solid and long-dashed lines) and inclination 
angle l (short-dashed and dotted lines) for model 9. The solid and short-dashed 
lines correspond to the outer disc section with r > 0.5, and the long-dashed and 
dotted lines correspond to the inner disc section with r < 0.3. 

Figure 9: The precession angle j3p (solid and long-daslicd lines) and inclination 
angle i (short-dashed and dotted lines) for model 2. The solid and short-dashed 
lines correspond to the outer disc section with r > 0.5, and the long-dashed and 
dotted lines correspond to the inner disc section with r < 0.3. 

Figure 10: Projection plots in orbital plane coordinates for model 4 at time 
t = 105.1. 

Figure 11: The precession angle /3p (dashed line) and inclination angle t (solid 
line) for model 4. 

Figure 12: Projection plots in orbital plane coordinates for model 11 at time 
t = 147.3. 
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Figure 13: Projection plots in orbital plane coordinates for model 6 at time 
t = 84.0. 



Figure 14: Projection plots in orbital plane coordinates for model 7 at time 
t = 328.5. 



Figure 15: Projection plots in orbital plane coordinates for model 8 at time 
t = 228.0. 



